C
C--ANDUCT SOLVES THE FREESTREAM VELOCITY GRADIENT EQUATION FOR AN
C--ANNULAR DUCT WITH VARIABLE WHIRL AND STAGNATION PRESSURE.
C
      COMMON NREAD,NWRIT
      DIMENSION STRFN (50),QDIST(50),ZLAMDA(50),TIP(50),RHOIP(50)
      DIMENSION R(101),ZLMINT(101),TINT(101),RHOINT(101),SAMP(101),     -
     1   A1(101),A2(101),A3(101),E(101),F(101),G1(101),DELPRS(101),     -
     2   DELLAM(101),DELTIP(101),RCARDQ(101),WSUBM(101),U(101),Q(101),  -
     3   AL(101),CURV(101),B1(101),W(101),BETA(101),PRES(101),UNEW(101),-
     4   WWCR(101),AAA(101),BBB(101)
      DIMENSION WHUB(100),WTIP(100),WWCRH(100),WWCRT(100),WMH(100),     -
     1   WMT(100),BETH(100),BETT(100),PRSH(100),PRST(100),HDIST(100),   -
     2   TDIST(100)
      DIMENSION TITLE(20),STITLE(20)
      DATA SBAND/.01/
      NREAD = 5
      NWRIT = 6
      PI = 3.1415927
      RTOLER = 1.E-3
    5 ISTAT = 0
      HDIST(1) = 0.
      TDIST(1) = 0.
      READ (NREAD,1000) TITLE
   10 WRITE(NWRIT,1100)
      ISTAT = ISTAT+1
      WRITE(NWRIT,1110) TITLE
      READ (NREAD,1000) STITLE
      WRITE(NWRIT,1110) STITLE
      READ (NREAD,1010) NHT,LSFR,IPRINT,NEXT
      WRITE(NWRIT,1120) NHT,LSFR,IPRINT,NEXT
      READ (NREAD,1020) GAM,AR,ZMSFL
      WRITE(NWRIT,1130) GAM,AR,ZMSFL
      READ (NREAD,1020) RHUB,RTIP,ZHUB,ZTIP,CURVH,CURVT,ALH,ALT
      WRITE(NWRIT,1140) RHUB,RTIP,ZHUB,ZTIP,CURVH,CURVT,ALH,ALT
      IF(LSFR.NE.0) GO TO 20
      READ (NREAD,1020) (STRFN(I),I=1,NHT)
      WRITE(NWRIT,1150) (STRFN(I),I=1,NHT)
      GO TO 30
   20 READ (NREAD,1020) (QDIST(I),I=1,NHT)
      WRITE(NWRIT,1160) (QDIST(I),I=1,NHT)
   30 READ (NREAD,1020) (ZLAMDA(I),I=1,NHT)
      WRITE(NWRIT,1170) (ZLAMDA(I),I=1,NHT)
      READ (NREAD,1020) (TIP(I),I=1,NHT)
      WRITE(NWRIT,1180) (TIP(I),I=1,NHT)
      READ (NREAD,1020) (RHOIP(I),I=1,NHT)
      WRITE(NWRIT,1190) (RHOIP(I),I=1,NHT)
C
C--ALL INPUT HAS BEEN READ
C
C
C--CALCULATE CONSTANTS AND INTERPOLATE AS NECESSARY, CONVERT ALH + ALT
C     TO RADIANS
C
      CP = AR/(GAM-1.)*GAM
      EXPON = 1./(GAM-1.)
      TGRGP1 = 2.*GAM*AR/(GAM+1.)
      PSI = ATAN2(ZHUB-ZTIP,RTIP-RHUB)
      QHT = SQRT((RTIP-RHUB)**2+(ZTIP-ZHUB)**2)
      DELQ = QHT/100.
      ALH = ALH/180.*PI
      ALT = ALT/180.*PI
      DELAL = (ALT-ALH)/100.
      JZ = 1
C
C--INITIALIZE AND MAKE INITIAL GUESS FOR WMHUB
C
      NCOUNT = 0
      NADD = 0
      NSUB = 0
      ZMXFLO = 0.
      ZMNFLO = 1.E10
      WMHUB = ZMSFL/RHOIP(1)/PI/(RHUB+RTIP)/QHT
      WCR = SQRT(TGRGP1*TIP(1))
      WMHUB = AMIN1(WMHUB,WCR)
      DELMAX = WMHUB/20.
      WMMIN = -1.E50
      DO 40 I=1,101
      R(I) = RHUB+FLOAT(I-1)/100.*(RTIP-RHUB)
      Q(I) = FLOAT(I-1)*DELQ
      U(I) = FLOAT(I-1)/100.*(R(I)+RHUB)/(RTIP+RHUB)
      AL(I) = ALH+FLOAT(I-1)*DELAL
   40 CURV(I) = CURVH+FLOAT(I-1)/100.*(CURVT-CURVH)
      IF(LSFR.NE.0) GO TO 60
C
C--STREAM FUNCTION GIVEN AS INPUT(LSFR = 0)
C
   50 CALL SPLINT(STRFN,ZLAMDA,NHT,U,101,ZLMINT,AAA,BBB)
      CALL SPLINT(STRFN,TIP   ,NHT,U,101,TINT  ,AAA,BBB)
      CALL SPLINT(STRFN,RHOIP ,NHT,U,101,RHOINT,AAA,BBB)
      GO TO 70
C
C--RADIUS GIVEN AS INPUT(LSFR = 1)
C
   60 CALL SPLINT(QDIST,ZLAMDA,NHT,Q,101,ZLMINT,AAA,BBB)
      CALL SPLINT(QDIST,TIP,NHT,Q,101,TINT,AAA,BBB)
      CALL SPLINT(QDIST,RHOIP,NHT,Q,101,RHOINT,AAA,BBB)
C
C--CALCULATE QUANTITIES WHICH ARE INDEPENDENT OF W-SUB-M
C
   70 PRSINT = RHOINT(1)*AR*TINT(1)
      DO 80 I=1,101
      SAMP(I) = SIN(AL(I)-PSI)
      CAMP = COS(AL(I)-PSI)
      A1(I) = CAMP*CURV(I)
      A2(I) = SAMP(I)/CAMP*CURV(I)
      A3(I) = -SIN(AL(I))/R(I)
      B1(I) = -SAMP(I)/CAMP
      E(I) = -ZLMINT(I)/R(I)/R(I)
      F(I) = (ZLMINT(I)/R(I))**2/2./TINT(I)
      G1(I) = AR/PRSINT
      RCARDQ(I) = RHOINT(I)*CAMP*2.*PI*R(I)*DELQ
      IF(I.GE.101) GO TO 80
      PRSNXT = RHOINT(I+1)*AR*TINT(I+1)
      DELPRS(I) = PRSNXT-PRSINT
      PRSINT = PRSNXT
      DELLAM(I) = ZLMINT(I+1)-ZLMINT(I)
      DELTIP(I) = TINT(I+1)-TINT(I)
   80 CONTINUE
C
C--SOLVE VELOCITY GRADIENT EQUATION ON GIVEN STRAIGHT LINE FROM HUB
C     TO TIP
C
C--START OR RESTART ITERATION PROCEDURE WITH IND = 1
  100 IND = 1
C--STATEMENTS 110 TO 140 PERFORM A SINGLE ITERATION
  110 WSUBM(1) = WMHUB
      NCOUNT = NCOUNT+1
      MSUP = 0
      MSON = 0
      ZMSFLE = 0.
      WTH2= (ZLMINT(1)/R(1))**2
      WM2 = WSUBM(1)**2
      TEMP = TINT(1)-(WTH2+WM2)/2./CP
      IF(TEMP.LE.0..AND.WSUBM(1).LT.0.) GO TO 150
      IF(TEMP.LE.0.) GO TO 160
      TTIP = TEMP/TINT(1)
      RVA = TTIP**EXPON*WSUBM(1)*RCARDQ(1)
      DO 140 I=1,100
      GRT = GAM*AR*TEMP
      IF(WM2.GT.GRT) MSUP = 1
      DENOM = 1.-WM2/GRT
      IF(ABS(DENOM).GT.SBAND) GO TO 120
      MSON = 1
      DENOM = SBAND
  120 ACOEF = A1(I)+SAMP(I)/DENOM*(A2(I)+A3(I)*(1.+WTH2/GRT))
      WAS = WSUBM(I)*(1.+ACOEF*DELQ+B1(I)/DENOM*DELAL+DELTIP(I)/2./     -
     1   TINT(I))
      WAS = WAS+(E(I)*DELLAM(I)+F(I)*DELTIP(I)+G1(I)*TEMP*DELPRS(I))/   -
     1   WSUBM(I)
      IF(WAS.LT.ABS(WMHUB/1000.)) GO TO 150
      WM2 = WAS**2
      WTH2 = (ZLMINT(I+1)/R(I+1))**2
      TEMP = TINT(I+1)-(WTH2+WM2)/2./CP
      IF(TEMP.LE.0..AND.WAS.LT.0.) GO TO 150
      IF(TEMP.LE.0.) GO TO 160
      GRT = GAM*AR*TEMP
      IF(WM2.GT.GRT) MSUP = 1
      DENOM = 1.-WM2/GRT
      IF(ABS(DENOM).GT.SBAND) GO TO 130
      MSON = 1
      DENOM = SBAND
  130 ACOEF = A1(I+1)+SAMP(I+1)/DENOM*(A2(I+1)+A3(I+1)*(1.+WTH2/GRT))
      WASS = WSUBM(I)+WAS*(ACOEF*DELQ+B1(I+1)/DENOM*DELAL+DELTIP(I)/2./ -
     1   TINT(I+1))
      WASS = WASS+(E(I+1)*DELLAM(I)+F(I+1)*DELTIP(I)+G1(I+1)*TEMP*      -
     1   DELPRS(I))/WAS
      IF(WASS.LT.ABS(WMHUB/1000.)) GO TO 150
      WSUBM(I+1) = (WAS+WASS)/2.
      WM2 = WSUBM(I+1)**2
      TEMP = TINT(I+1)-(WTH2+WM2)/2./CP
      IF(TEMP.LE.0..AND.WSUBM(I+1).LT.0.) GO TO 150
      IF(TEMP.LE.0.) GO TO 160
      IF(WSUBM(I+1).LT.0.) GO TO 150
      TTIP = TEMP/TINT(I+1)
      RVAS = TTIP**EXPON*WSUBM(I+1)*RCARDQ(I+1)
      ZMSFLE = (RVA+RVAS)/2.+ZMSFLE
      UNEW(I+1) = ZMSFLE/ZMSFL
  140 RVA = RVAS
      IF(IND.GE.6.AND.ABS(ZMSFL-ZMSFLE).LE.ZMSFL*RTOLER) GO TO 200
      ZMXFLO = AMAX1(ZMSFLE,ZMXFLO)
      ZMNFLO = AMIN1(ZMSFLE,ZMNFLO)
      IF(WMMIN.LT.-1.E40) WMMAX = WMHUB
      IF(WMMIN.LT.-1.E40) WMMIN = WMHUB
      WMMAX = AMAX1(WMHUB,WMMAX)
      WMMIN = AMIN1(WMHUB,WMMIN)
      CALL CONTIN(WMHUB,ZMSFLE,IND,JZ,ZMSFL,DELMAX)
      IF(WMHUB.EQ.0.) GO TO 300
      IF(IND.LT.10) GO TO 110
C--IND = 10 INDICATES CHOKED FLOW
      IF(IND.EQ.10) GO TO 200
C--IND = 11 INDICATES NO SOLUTION FOUND IN 100 ITERATIONS
      GO TO 300
  150 WMHUB = WMHUB+.501*DELMAX
      NADD = NADD+1
      IF(NCOUNT.LT.1000) GO TO 100
      GO TO 300
  160 WMHUB = WMHUB-.499*DELMAX
      NSUB = NSUB+1
      IF(NCOUNT.LT.1000) GO TO 100
      GO TO 300
C
C--CONVERGED INNER ITERATION
C
  200 CONTINUE
      IF(LSFR.NE.0) GO TO 220
      ERRMAX = 0.
      DO 210 I=2,101
      IF(LSFR.EQ.0) ERRMAX = AMAX1(ERRMAX,ABS(UNEW(I)-U(I)))
  210 U(I) = UNEW(I)
      IF(ERRMAX.GT.RTOLER.AND.NCOUNT.LT.1000) GO TO 50
C
C--CONVERGED OUTER ITERATION - PRINT DETAILED OUTPUT IF REQUESTED
C
  220 IF(IPRINT.NE.0) WRITE(NWRIT,1100)
      IF(IPRINT.NE.0) WRITE(NWRIT,1110) STITLE
      IF(ERRMAX.GT.RTOLER) WRITE(NWRIT,2030) ERRMAX
      IF(IND.EQ.10) WRITE(NWRIT,2000) ZMSFLE
      IF(MSUP.EQ.1) WRITE(NWRIT,2010)
      IF(MSON.EQ.1) WRITE(NWRIT,2020)
      DO 230 I=1,101
      WTH = ZLMINT(I)/R(I)
      WSQ = WSUBM(I)**2+WTH**2
      W(I) = SQRT(WSQ)
      WWCR(I) = W(I)/SQRT(TGRGP1*TINT(I))
      BETA(I) = ATAN2(WTH,WSUBM(I))/PI*180.
      PRSINT = RHOINT(I)*AR*TINT(I)
      TTIP = 1.-WSQ/2./CP/TINT(I)
  230 PRES(I) = PRSINT*TTIP**(GAM*EXPON)
      IF(ISTAT.GT.1) HDIST(ISTAT) = HDIST(ISTAT-1)+SQRT((RHUB-RHLAST)**2-
     1   +(ZHUB-ZHLAST)**2)
      IF(ISTAT.GT.1) TDIST(ISTAT) = TDIST(ISTAT-1)+SQRT((RTIP-RTLAST)**2-
     1   +(ZTIP-ZTLAST)**2)
      RHLAST = RHUB
      RTLAST = RTIP
      ZHLAST = ZHUB
      ZTLAST = ZTIP
      WHUB(ISTAT) = W(1)
      WTIP(ISTAT) = W(101)
      WWCRH(ISTAT) = WWCR(1)
      WWCRT(ISTAT) = WWCR(101)
      WMH(ISTAT) = WSUBM(1)
      WMT(ISTAT) = WSUBM(101)
      BETH(ISTAT) = BETA(1)
      BETT(ISTAT) = BETA(101)
      PRSH(ISTAT) = PRES(1)
      PRST(ISTAT) = PRES(101)
      IF(IPRINT.EQ.0) GO TO 500
      WRITE(NWRIT,1200) (I,W(I),WWCR(I),WSUBM(I),BETA(I),PRES(I),U(I)   -
     1   ,I=1,101,2)
      GO TO 500
C
C--PRINT ERROR MESSAGES WITH INFORMATION IF A SATISFACTORY SOLUTION
C     CANNOT BE OBTAINED
  300 CONTINUE
      WRITE(NWRIT,1100)
      WRITE(NWRIT,1110) STITLE
      IF(IND.EQ.11) WRITE(NWRIT,2040)
      IF(NCOUNT.GE.1000) WRITE(NWRIT,2050)
      WRITE(NWRIT,2060) ZMXFLO,ZMNFLO,WMMAX,WMMIN,NCOUNT
      IF(NSUB.GT.0) WRITE(NWRIT,2070) NSUB
      IF(NADD.GT.0) WRITE(NWRIT,2080) NADD
  500 IF(NEXT.EQ.1.AND.ISTAT.LT.100) GO TO 10
C
C--PRINT SUMMARY FOR A GIVEN GEOMETRY
C
      IF(NEXT.EQ.1.AND.ISTAT.GE.100) WRITE(NWRIT,2110)
      WRITE(NWRIT,1100)
      WRITE(NWRIT,1110) TITLE
      WRITE(NWRIT,2130) ISTAT
      WRITE(NWRIT,2120) (HDIST(I),WHUB(I),WWCRH(I),WMH(I),BETH(I),      -
     1   PRSH(I),I=1,ISTAT)
      WRITE(NWRIT,2140) ISTAT
      WRITE(NWRIT,2120) (TDIST(I),WTIP(I),WWCRT(I),WMT(I),BETT(I),      -
     1   PRST(I),I=1,ISTAT)
C
C--CONTINUE WITH NEXT CASE IF ANY
C
      IF(NEXT.EQ.0) STOP
      IF(NEXT.EQ.-1) GO TO 5
      ISTAT = 0
      GO TO 10
C--FORMAT STATEMENTS
 1000 FORMAT(20A4)
 1010 FORMAT(16I5)
 1020 FORMAT(8F10.0)
 1100 FORMAT(1H1)
 1110 FORMAT(20X,20A4)
 1120 FORMAT(1HL,5X,3HNHT,3X,4HLSFR,1X,6HIPRINT,3X,4HNEXT/2X,16I7)
 1130 FORMAT(5X,3HGAM,13X,2HAR,12X,5HZMSFL/1X,8G16.7)
 1140 FORMAT(1HL,4X,4HRHUB,12X,4HRTIP,12X,4HZHUB,12X,4HZTIP,11X,5HCURVH,-
     1   11X,5HCURVT,12X,3HALH,13X,3HALT/1X,8G16.7)
 1150 FORMAT(1HL,10X,11HSTRFN ARRAY/(1X,8G16.7))
 1160 FORMAT(1HL,10X,11HQDIST ARRAY/(1X,8G16.7))
 1170 FORMAT(1HL,10X,12HZLAMDA ARRAY/(1HX,8G16.7))
 1180 FORMAT(1HL,10X,9HTIP ARRAY/(1X,8G16.7))
 1190 FORMAT(1HL,10X,11HRHOIP ARRAY/(1X,8G16.7))
 1200 FORMAT(1HK,2X,1HI,7X,1HV,13X,5HV/VCR,11X,5HVSUBM,12X,4HBETA,6X,   -
     1   15HSTATIC PRESSURE,3X,15HSTREAM FUNCTION/(2X,I3,6G16.7))
 2000 FORMAT(1HL,120(1H*)/1HL,24X,38HTHE PASSAGE IS CHOKED AT THIS STATI-
     1ON /26X,24HTHE CHOKING MASS FLOW IS ,G12.4/1HL,120(1H*))
 2010 FORMAT(1HL,10X,56HSUPERSONIC MERIDIONAL VELOCITY COMPONENT AT THIS-
     1 STATION)
 2020 FORMAT(1HL,120(1H*)/1HL,30X,58HSONIC MERIDIONAL VELOCITY COMPONENT-
     1 OCCURS AT THIS STATION/40X,41HTHIS MAY RESULT IN AN INACCURATE SO-
     2LUTION/1HL,120(1H*))
 2030 FORMAT(1HL,120(1H*)/1HL,20X,83HA FULLY CONVERGED SOLUTION COULD NO-
     1T BE OBTAINED IN 1000 ITERATIONS AT THIS STATION/30X,30HTHE STREAM-
     2 FUNCTION CHANGED BY,F6.3,32H BETWEEN THE LAST TWO ITERATIONS/1HL,-
     3   120(1H*))
 2040 FORMAT(10X,44HNO SOLUTION COULD BE FOUND IN 100 ITERATIONS)
 2050 FORMAT(10X,69HITERATION PROCEDURE HAD TO BE RESTARTED TO AVOID NEG-
     1ATIVE TEMPERATURE/15X,67HRESTART PROCEDURE WAS ABORTED AFTER 1000 -
     2TOTAL NUMBER OF ITERATIONS)
 2060 FORMAT(1HL,64HTHE MAXIMUM MASS FLOW FOR WHICH A SOLUTION COULD BE -
     1OBTAINED WAS,G16.7/ 1X,64HTHE MINIMUM MASS FLOW FOR WHICH A SOLUTI-
     2ON COULD BE OBTAINED WAS,G16.7/1HL,80HTHE MAXIMUM VALUE OF VSUBM A-
     3T THE HUB FOR WHICH A SOLUTION COULD BE OBTAINED WAS,G16.7/1X,80HT-
     4HE MINIMUM VALUE OF VSUBM AT THE HUB FOR WHICH A SOLUTION COULD BE-
     5 OBTAINED WAS,G16.7/1HL,34HTHE TOTAL NUMBER OF ITERATIONS WAS,I7)
 2070 FORMAT(1HL,6HNSUB =,I7)
 2080 FORMAT(1HL,6HNADD =,I7)
 2110 FORMAT(1HL,120(1H*)/1HL,20X,52HTHE LIMIT OF 100 STATIONS PER CASE -
     1HAS BEEN EXCEEDED/25X,47HOUTPUT IS GIVEN ONLY FOR THE FIRST 100 ST-
     2ATIONS/1HL,120(1H*))
 2120 FORMAT(1HL,5X,8HDISTANCE,12X,1HV,11X,5HV/VCR,11X,7HV-SUB-M,11X,   -
     1   4HBETA,7X,15HSTATIC PRESSURE/(2X,6G16.6))
 2130 FORMAT(1HL,19HHUB OUTPUT DATA FOR,I4,9H STATIONS)
 2140 FORMAT(1HL,22HSHROUD OUTPUT DATA FOR,I4,9H STATIONS)
      END
      SUBROUTINE CONTIN(XEST,YCALC,IND,JZ,YGIV,XDEL)
C
C--CONTIN CALCULATES AN ESTIMATE OF THE RELATIVE FLOW VELOCITY
C--FOR USE IN THE VELOCITY GRADIENT EQUATION
C
      DIMENSION X(3),Y(3)
      NCALL = NCALL+1
      IF (IND.NE.1.AND.NCALL.GT.100) GO TO 160
      GO TO (10,30,40,50,60,110,150),IND
C--FIRST CALL
   10 NCALL = 1
      XORIG = XEST
      IF (YCALC.GT.YGIV.AND.JZ.EQ.1) GO TO 20
      IND = 2
      Y(1) = YCALC
      X(1) = 0.
      XEST = XEST+XDEL
      RETURN
   20 IND = 3
      Y(3) = YCALC
      X(3) = 0.
      XEST = XEST-XDEL
      RETURN
C--SECOND CALL
   30 IND = 4
      Y(2) = YCALC
      X(2) = XEST-XORIG
      XEST = XEST+XDEL
      RETURN
   40 IND = 5
      Y(2) = YCALC
      X(2) = XEST-XORIG
      XEST = XEST-XDEL
      RETURN
C--THIRD OR LATER CALL - FIND SUBSONIC OR SUPERSONIC SOLUTION
   50 Y(3) = YCALC
      X(3) = XEST-XORIG
      GO TO 70
   60 Y(1) = YCALC
      X(1) = XEST-XORIG
   70 IF (YGIV.LT.AMIN1(Y(1),Y(2),Y(3))) GO TO (120,130),JZ
   80 IND = 6
      CALL PABC(X,Y,APA,BPB,CPC)
      DISCR = BPB**2-4.*APA*(CPC-YGIV)
      IF (DISCR.LT.0.) GO TO 140
      IF (ABS(400.*APA*(CPC-YGIV)).LE.BPB**2) GO TO 90
      XEST = -BPB-SIGN(SQRT(DISCR),APA)
      IF (JZ.EQ.1.AND.APA.GT.0..AND.Y(3).GT.Y(1)) XEST = -BPB+          -
     1SQRT(DISCR)
      IF (JZ.EQ.2.AND.APA.LT.0.) XEST = -BPB-SQRT(DISCR)
      XEST = XEST/2./APA
      GO TO 100
   90 IF (JZ.EQ.2.AND.BPB.GT.0.) GO TO 130
      ACB2 = APA/BPB*(CPC-YGIV)/BPB
      IF (ABS(ACB2).LE.1.E-8) ACB2=0.
      XEST = -(CPC-YGIV)/BPB*(1.+ACB2+2.*ACB2**2)
  100 IF (XEST.GT.X(3)) GO TO 130
      IF (XEST.LT.X(1)) GO TO 120
      XEST = XEST+XORIG
      RETURN
C--FOURTH OR LATER CALL - NOT CHOKED
  110 IF(XEST-XORIG.GT.X(3)) GO TO 130
      IF(XEST-XORIG.LT.X(1)) GO TO 120
      Y(2) = YCALC
      X(2) = XEST-XORIG
      GO TO 70
C--THIRD OR LATER CALL - SOLUTION EXISTS,
C--BUT RIGHT OR LEFT SHIFT REQUIRED
  120 IND = 5
C--LEFT SHIFT
      XEST = X(1)-XDEL+XORIG
      XOSHFT = XEST-XORIG
      XORIG = XEST
      Y(3) = Y(2)
      X(3) = X(2)-XOSHFT
      Y(2) = Y(1)
      X(2) = X(1)-XOSHFT
      RETURN
  130 IND = 4
C--RIGHT SHIFT
      XEST = X(3)+XDEL+XORIG
      XOSHFT = XEST-XORIG
      XORIG = XEST
      Y(1) = Y(2)
      X(1) = X(2)-XOSHFT
      Y(2) = Y(3)
      X(2) = X(3)-XOSHFT
      RETURN
C--THIRD OR LATER CALL - APPEARS TO BE CHOKED
  140 XEST = -BPB/2./APA
      IND = 7
      IF (XEST.LT.X(1)) GO TO 120
      IF(XEST.GT.X(3)) GO TO 130
      XEST = XEST+XORIG
      RETURN
C--FOURTH OR LATER CALL - PROBABLY CHOKED
  150 IF (YCALC.GE.YGIV) GO TO 110
      IND = 10
      RETURN
C--NO SOLUTION FOUND IN 100 ITERATIONS
  160 IND = 11
      RETURN
      END
      SUBROUTINE PABC(X,Y,A,B,C)
C
C--PABC CALCULATES COEFFICIENTS A,B,C OF THE PARABOLA
C--Y=A*X**2+B*X+C, PASSING THROUGH THE GIVEN X,Y POINTS
C
      DIMENSION X(3),Y(3)
      C1 = X(3)-X(1)
      C2 = (Y(2)-Y(1))/(X(2)-X(1))
      A = (C1*C2-Y(3)+Y(1))/C1/(X(2)-X(3))
      B = C2-(X(1)+X(2))*A
      C = Y(1)-X(1)*B-X(1)**2*A
      RETURN
      END
      SUBROUTINE SPLINT (X,Y,N,Z,MAX,YINT,DYDX,D2YDX2)
C
C--SPLINT CALCULATES INTERPOLATED POINTS AND DERIVATIVES
C--FOR A SPLINE CURVE
C--END CONDITION - SECOND DERIVATIVES AT END POINTS ARE
C--SDR1 AND SDRN TIMES SECOND DERIVATIVES AT ADJACENT POINTS
C
      COMMON NREAD,NWRIT
      DIMENSION X(N),Y(N),Z(MAX),YINT(MAX),DYDX(MAX),D2YDX2(MAX)
      DIMENSION G(101),SB(101),EM(101)
      IERR = 0
      SDR1 = .5
      SDRN = .5
      TOLER= ABS(X(N)-X(1))/FLOAT(N)*1.E-5
      C = X(2)-X(1)
      IF (C.EQ.0.) GO TO 130
      SB(1) = -SDR1
      G(1) = 0.
      NO = N-1
      IF (NO.LE.0) GO TO 140
      IF (NO.EQ.1) GO TO 20
      DO 10 I=2,NO
      A = C
      C = X(I+1)-X(I)
      IF (A*C.EQ.0.) GO TO 130
      IF (A*C.LT.0.) IERR = 1
      W = 2.*(A+C)-A*SB(I-1)
      SB(I) = C/W
      F = (Y(I+1)-Y(I))/C-(Y(I)-Y(I-1))/A
   10 G(I) = (6.*F-A*G(I-1))/W
   20 EM(N) = SDRN*G(N-1)/(1.+SDRN*SB(N-1))
      DO 30 I=2,N
      K = N+1-I
   30 EM(K) = G(K)-SB(K)*EM(K+1)
      IF (MAX.LE.0) RETURN
C
      ENTRY SPLENT (Z,MAX,YINT,DYDX,D2YDX2)
      DO 120 I=1,MAX
      K=2
      IF (ABS(Z(I)-X(1)).LT.TOLER) GO TO 40
      IF (Z(I).GT.2.0*X(1)-X(2)) GO TO 50
      GO TO 80
   40 YINT(I) = Y(1)
      SK = X(K)-X(K-1)
      GO TO 110
   50 IF (ABS(Z(I)-X(K)).LT.TOLER) GO TO 60
      IF (Z(I).GT.X(K)) GO TO 70
      GO TO 100
   60 YINT(I) = Y(K)
      SK = X(K)-X(K-1)
      GO TO 110
   70 IF (K.GE.N) GO TO 90
      K = K+1
      GO TO 50
   80 S2 = X(2)-X(1)
      Y0 = EM(1)*S2**2+2.*Y(1)-Y(2)
      DYDX(I) = (Y(2)-Y(1))/S2-7.*EM(1)/6.*S2
      YINT(I) = Y0+DYDX(I)*(Z(I)-X(1)+S2)
      D2YDX2(I) = 0.
      GO TO 120
   90 IF (Z(I).LT.2.*X(N)-X(N-1)) GO TO 100
      SN = X(N)-X(N-1)
      YNP1 = EM(N)*SN**2+2.*Y(N)-Y(N-1)
      DYDX(I) = (Y(N)-Y(N-1))/SN+7.*EM(N)/6.*SN
      YINT(I) = YNP1+DYDX(I)*(Z(I)-X(N)-SN)
      D2YDX2(I) = 0.
      GO TO 120
  100 SK = X(K)-X(K-1)
      YINT(I) = EM(K-1)*(X(K)-Z(I))**3/6./SK  +EM(K)*(Z(I)-X(K-1))**3/6.-
     1  /SK+(Y(K)/SK  -EM(K)*SK  /6.)*(Z(I)-X(K-1))+(Y(K-1)/SK  -EM(K-1)-
     2  *SK/6.)*(X(K)-Z(I))
  110 DYDX(I)=-EM(K-1)*(X(K)-Z(I))**2/2.0/SK  +EM(K)*(X(K-1)-Z(I))**2/2.-
     1   /SK+(Y(K)-Y(K-1))/SK  -(EM(K)-EM(K-1))*SK/6.
      D2YDX2(I) = EM(K)-(X(K)-Z(I))/SK*(EM(K)-EM(K-1))
  120 CONTINUE
      IF (IERR.EQ.0) RETURN
  130 WRITE(NWRIT,1000)
      WRITE(NWRIT,1020) N,(X(I),Y(I),I=1,N)
      IF (IERR.EQ.0) STOP
      IERR = 0
      RETURN
  140 WRITE(NWRIT,1010)
      WRITE(NWRIT,1020) N,(X(I),Y(I),I=1,N)
      STOP
 1000 FORMAT (1H1,10X,44HSPLINT ERROR -- ONE OF THREE POSSIBLE CAUSES/  -
     117X,51H1.  ADJACENT X POINTS ARE DUPLICATES OF EACH OTHER./       -
     217X,38H2.  SOME X POINTS ARE OUT OF SEQUENCE./                    -
     317X,32H3.  SOME X POINTS ARE UNDEFINED.)
 1010 FORMAT (1H1,10X,62HSPLINT ERROR -- NUMBER OF SPLINE POINTS GIVEN I-
     1S LESS THAN TWO)
 1020 FORMAT (//17X,18HNUMBER OF POINTS =,I4//17X,8HX  ARRAY,6X,8HY  ARR-
     1AY/(17X,2G13.5))
      END
